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Abstract 

The rapid progress of genomic technologies has been providing new opportunities to address the need of maturity-onset 
diabetes of the young (IVIODY) molecular diagnosis. However, whether a new mutation causes IVIODY can be questionable. 
A number of In silico methods have been developed to predict functional effects of rare human mutations. The purpose of 
this study is to compare the performance of different bioinformatics methods in the functional prediction of 
nonsynonymous mutations in each IVIODY gene, and provides reference matrices to assist the molecular diagnosis of 
MODY. Our study showed that the prediction scores by different methods of the diabetes mutations were highly correlated, 
but were more complimentary than replacement to each other. The available in silico methods for the prediction of diabetes 
mutations had varied performances across different genes. Applying gene-specific thresholds defined by this study may be 
able to increase the performance of in silico prediction of disease-causing mutations. 
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introduction 

To date, a number of methods have been developed to predict 
fiinctional effects of rare human mutations based on the impact of 
protein function and/or evolutionary conservation [1-3]. These 
methods are valuable to assist the diagnosis of monogenic 
inheritance diseases. In the area of diabetes, there is a common 
monogenic form, i.e. maturity-onset diabetes of the young 
(MODY). MODY accounts for -1% to 5% of all cases of 
diabetes, while it is mainly seen in young adults (£25 years old) 
[4] . As an autosomal dominant inherited form of diabetes, MODY 
is caused by gene mutations leading to insufficient insulin 
production without or with minimal insulin resistance [5]. To 
date, at least 13 genes have been identified with mutations that 
cause MODY, i.e. HNF4A (MODYl) [6], GCK (MODY2) [7], 
HNFIA (MODY3) [6], PDXl (MODY4) [8], HNFIB (MODY5) 
[9], NEURODl (MODY6) [10], KLFll (MODY7) [11], CEL 
(MODY8) [12], PAX4 (MODY9) [13], INS (MODY 10) [14], 
5L^(MODYll) [15],^fiCC<§(MODY12) [16,17], andKCNJll 
(MODYl 3) [16,17]. 

MODY caused by different gene mutations may have different 
severities of diabetes and different drug responses to diabetes 
medications [18]. For example, MODY2 (accounts for ~20% of 
all MODY cases [19]) caused by GCK gene mutations tends to 
have mild hyperglycaemia without obvious glycosuria. Patients 
with MODY2 are often asymptomatic [20] or only identified in 



women during pregnancy and diagnosed as gestational diabetes 
[21]. Most patients with MODY2 can have blood glucose 
satisfactorily controlled by diet therapy and don't need hypogly- 
cemic medication [22,23]. In contrast, MODY3 caused by 
HNFIA mutations, the most common type of MODY that 
accounts for ~63% of all MODY cases [19], tends to have obvious 
glycosuria because of impaired glucose-stimulated insulin secretion 
[24], as well as decreased renal threshold for glucose [25]. 
MODY3 patients tend to have good response to sulphonylurea 
treatment and don't rely on insulin therapy [18,26]. Because of the 
implications of pharcogenetics and personalized medicine, molec- 
ular diagnosis of MODY has clinical importance for clinical 
decision and for genetic counseling [18,26]. However, because of 
unavailability and expense of MODY molecular diagnosis, it is not 
vmcommon that MODY patients are classified as type 2 diabetes 
[27,28] and occasionally as type 1 [29]. 

The rapid progress of advanced genomic technologies has been 
providing new opportunities to address the need of MODY 
molecular diagnosis. The identification of mutations in MODY 
genes by sequencing technologies will enable the molecular 
diagnosis of MODY, whereas a new issue is emerging. Most 
mutations causing MODY are nonsynonymous single-nucleotide 
mutations causing the change of an amino acid residue (according 
to The Human Gene Mutation Database [30], http:/ /www.hgmd. 
org/). High throughput sequencing technologies enable screening 
of a large number of patients and parallel sequencing of a large 
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Table 1. MODY or neonatal diabetes genes and mutations (n 


= 1091). 






Gene symbol 


Diabetic mutations* (n) 


Control mutations (n) 


Chromosome position 


Gene name 


HNF4A 


86 


57 


20ql2-ql3.1 


hepatocyte nuclear factor 4, alpha 


GCK 


479 


22 


7pl5.3-pl5.1 


glucokinase (hexokinase 4, maturity onset 
diabetes of the young 2) 


HNFIA 


324 


78 


12q24.2 


hepatocyte nuclear factor-1 (HNFl) homeobox A 


HNFIS 


36 


52 


17cen-q21.3 


hepatocyte nuclear factor-1 (HNFl) homeobox B 


INS 


41 


14 


llpl5.5 


insulin 


ABCC8 


64 


185 


llpl5.1 


ATP-binding cassette, sub-family C (CFTR/MRP), 
member 8 


KCNJU 


61 


65 


llpl5.1 


potassium inwardly-rectifying channel, 
subfamily J, member 11 



* Number of diabetes-causing single nucleotide mutations recorded in the Human Gene IWutation Database (HGMD) 2013.4 release (http://www.hgmd.org/) [30]. 
doi:l 0.1 371 /journal.pone.Ol 04452.t001 



number of genes. If a known MODY gene mutation is identified in sequencing technologies is likely to produce increased numbers 
a patient suspected of MODY, the molecular diagnosis of MODY of missense variants whose causative role in MODY can be 
can be established. However, the increased throughput of questionable. Biornformatics tools, e.g. SIFT (http:/ /sift.jcvi.org/) 



Table 2. Methods for function prediction for non-synonymous mutations*. 



Method 


Deleterious Threshold 


Algorithm 


PhyloP 


>1.6 


PhyloP calculates basewise conservation score from Multiz alignment [55] of 46 vertebrate species {ftp:// 
ccg.vital-it.ch/mga/hg 1 9/phylop/phylop.html) [38]. 


GERP++ RS 


>4.4 


GERP++ RS calculates site-specific "rejected substitutions" (RS) scores and to discover evolutionarily 
constrained elements based on maximum likelihood evolutionary rate estimation [39]. 


SiPhy 


>12.17 


SiPhy detects bases under selection from a multiple alignment data using a hidden Markov model, (http:// 
www.broadinstitute.org/genome_bio/siphy/) [40]. 


SIFT 


>0.95 


SIFT prediction is based on the degree of conservation of amino acid residues in sequence alignments 
derived from closely related sequences, collected through PSI-BLAST.(http://sift.jcvi.org/) [31]. 


PolyPhen-2 


>0.5 


PolyPhen-2 (Polymorphism Phenotyping v2) predicts the functional significance of an amino acid 
substitution by Naive Bayes classifier, using sequence-based and structure-based predictive features. HDIV, 
or HumDiv, identifies human damaging mutations by assuming differences between human proteins and 
their closely related mammalian homologs as non-damaging. HVAR, or HumVar, identifies human disease- 
causing mutations by assuming common human nsSNPs as non-damaging, (http://genetics.bwh. harvard. 
edu/pph2/index.shtml) [32]. 


LRT 


>0.999 


The likelihood ratio test (LRT) identifies conserved amino acid positions and deleterious mutations using a 
comparative genomics data set of multiple vertebrate species, (http://www.genetics.wustl.edu/jflab/ 
lrt_query.html) [41]. 


MutationTaster 


>0.5 


MutationTaster evaluates the disease-causing potential of DNA sequence alterations by Naive Bayes 
classifier, integrating information Of evolutionary conservation, splice-site changes, loss of protein features 
and changes that might affect the amount of mRNA from different biomedical databases and uses 
established analysis tools, (http://www.mutationtaster.org/) [42]. 


Mutation Assessor 


>0.65 


Mutation Assessor predicts the functional impact of amino-acid substitutions in proteins based on 
evolutionary conservation of the affected amino acid in protein homologs. (http://mutationassessor.org) 
[43]. 


PATH MM 


> 0.45 3 


The Functional Analysis through Hidden Markov Models (FATHMM) Predicts the functional consequences 
of cancer-associated amino acid substitutions using a model weighted for inherited disease mutations 
(http://fathmm.biocompute.org.uk) [44]. 


RadialSVM score 


>0.5 


RadialSVM score is an ensemble-based approach integrating multiple scoring systems (function prediction 
and conservation Score) by radial support vector machine (SVM) [3]. 


LR score 


>0.5 


LR score is an ensemble-based approach integrating multiple scoring systems (function prediction and 
conservation Score) by logistic regression (LR) [3]. 



* Extensively comparisons of these methods genome-widely have been studied byThusberg et al. [1,2] and Dong et al. [3]. Deleterious thresholds for PhyloP, GERP-F-F RS 

and SiPhy are defined according to the study by Dong et al. [3]. Deleterious thresholds for SIFT, LRT, MutationTaster, Mutation Assessor, FATHMM and RAdialSVM are 

based on converted scores used in dbNSFP version 2.3 [45,46] (converted score is designated as S^ and original score is designated as So): SIFT: Sc = I-Sq; LRT: 5^ = 1- So 

*0.5 if n<l, or Sc = So *0.5 if Cl> = 1; MutationTaster: Sc = So if the prediction is "A" or "D" or Sc = 1- So if the prediction is "N" or "P"; Mutation Assessor: 

Sc = (So -(-5.545))/(5.975-(-5.545)); FATHMM: Sc=l-(So -(- 16.1 3))/(10.64-(- 16.13)); RadialSVM: S^ = C\+ So/3.03993691 875303)*0.5 if predicted "D" and 

Sc = (l- So/-2.00575697514507)*0.5 otherwise. More details of the conversion can be found at http://dbnsfp.houstonbioinformatics.Org/dbNSFPzip/dbNSFP2.3.readme. 

txt. 

doi:l 0.1 371 /journal.pone.Ol 04452.t002 
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Table 4. Method comparisons for function prediction for non-synonymous mutations causing diabetes. 



Methods 


Missing Rate* 


False Negative Rate 


False Positive Rate** 


MCC 


PhyloP 


0% 


18% 


53% 


0.300 


GERP-H- RS 


0% 


21% 


52% 


0.281 


SiPhy 


0% 


16% 


51% 


0.342 


SIFT 


13% 


25% 


39% 


0.350 


PolyPhen-2 HDIV 


15% 


9% 


51% 


0.447 


PolyPhen-2 HVAR 


15% 


16% 


42% 


0.434 


LRT 


18% 


7% 


68% 


0.324 


MutationTaster 


3% 


3% 


77% 


0.333 


Mutation Assessor 


15% 


30% 


32% 


0.362 


FATHMM 


14% 


1% 


95% 


0.127 


RadialSVM score 


8% 


5% 


57% 


0.474 


LR score 


8% 


4% 


69% 


0.393 



* The missing rate refers to the percentage of mutations that a method is inapplicable; 

**The false positive rate was calculated by nonsynonymous single-nucleotide mutations in the diabetes genes acquired from the NHLBI GO Exome Sequencing Project 
(ESP) [33], the CHARGE Exome Sequencing Project [34,36], and the 1000 Genome Project [37], excluding mutations recorded in the HGIVID database. 
doi:l 0.1 371 /journal.pone.01 04452.t004 



[31] and PolyPhen (http://genetics.bwh.harvard.edu/pph2/ 
index.shtml) [32], are often used to assess the pathogenicity of a 
nonsynonymous mutation [27]. Because the limitations of m silico 
methods, the functional prediction of a nonsynonymous mutation 
lacks a gold standard. To date, a number of bioinformatics 
methods besides SIFT and PolyPhen, based on difiFerent 
algorithms, have been developed [1-3]. The purpose of this study 
is to compare the performances of different bioinformatics 
methods in the functional prediction of nonsynonymous mutations 
in each MODY gene, and to provide reference matrices to assist 
the molecular diagnosis of MODY. 

Methods 

Data source 

The diabetes mutation data analyzed in this study were 
acquired from the Human Gene Mutation Database (HGMD) 
2013.4 release (http://www.hgmd.org/) [30]. As the purpose of 
this study is to assess the prediction performances of different in 
silico methods for nonsynonymous single-nucleotide mutations, 
insertion/ deletion mutations (InDels) are not included in this 
study. Altogether, 1,130 nonsynonymous single-nucleotide muta- 
tions from 24 genes have been reported causing MODY or 
neonatal diabetes. Among these genes, 7 genes harbor more than 
30 single-nucleotide mutations within each gene with the total of 
1,091 diabetes mutations (Table 1), while the other 17 genes 
harbor £6 diabetes mutations in each gene. To enable statistical 
comparisons of different in silico methods across different genes, 
those 17 genes with £6 diabetes mutations were not involved in 
this study. Among the 1,091 mutations, 155 mutations from the 
genes ABCC8, GCK, INS, or KCNJll, have the phenotype of 
neonatal diabetes, either transient or permanent. The other 936 
mutations in the 7 genes have the phenotype of MODY. 

Control nonsynonymous single-nucleotide mutations in the 
diabetes genes were acquired from the NHLBI GO Exome 
Sequencing Project (ESP) [33,34], the ARIC samples [35] in the 
CHARGE Exome Sequencing Project [36], and the 1000 
Genome Project [37], excluding mutations recorded in the 
HGMD database. 



Functional prediction of nonsynonymous single- 
nucleotide mutations 

Eleven methods, including PhyloP [38], GERP+-F RS [39], 
SiPhy [40], SIFT [31], PolyPhen-2 [32], the Ukelihood ratio test 
(LRT) [41], MutationTaster [42], Mutation Assessor [43], 
FATHMM [44], RadialSVM score [3], and logistic regression 
(LR) score [3], were covered in the dbNSFP database [45,46] and 
compared in this study (Table 2). Among the 1,091 mutations 
involved in this analysis, 104 mutations from the genes GCK, 
HNFIA, HNFIB, HNF4A, and INS, are nonsense mutations, i.e. 
producing a pre-termination codon; two other mutations from the 
gene GCK replace a termination codon with an amino acid codon. 
For these mutations, the methods, PolyPhen-2 HDIV, PolyPhen-2 
HVAR, MutationAssessor, FATHMM, Radial SVM score, LR 
score are not apphcable to nonsense mutations or mutations 
assumed with highly damaging potential. Other methods, except 
MutationTaster, tend to have higher error rates (false negative 
rates, FNR), compared to the prediction of amino acid substitution 
mutations, i.e. SIFT FNR = 72%, GERP+-F RS FNR = 40%, 
PhyloP FNR =35%, SiPhy FNR = 26%, LRT FNR = 21%, and 
MutationTaster FNR = 3%, for the prediction of nonsense 
mutations. Compared with amino acid substitutions, the assess- 
ment of the functional effect of nonsense mutations tends to be less 
of an issue. The final analysis of this study involved 985 
nonsynonymous single-nucleotide mutations. The quantitative 
performances of these methods were compared by the Spearman's 
rank correlation test [47] and the ANOVA test using the IBM 
SPSS Statistics 19 software ((IBM SPSS Inc., Chicago, IE, USA), 
To re-define gene-specific thresholds of deleterious mutations, the 
receiver operator characteristic (ROC) analysis was calculated by 
the sensitivity and specificity values of screening series of cutoffs of 
each method for each gene. A redefined threshold was identified 
using the maximum Matthews correlation coefficient (MCC) [48]. 

Results and Discussion 

In our analysis, the prediction scores by different methods of the 
diabetes mutations are highly correlated (Table 3). The highest 
correlations are seen between RadialSVM score and LR score 
(r = 0.957), PolyPhen-2 HDIV and PolyPhen-2 HVAR (r = 0.89). 
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and phyloP and GERP++ RS (r = 0.871), while the other 
correlations have r<0.80. Therefore, in spite of the high statistical 
significance of correlations between different methods, different 
methods may not he able to replace each other except for the 
above three pairs. Especially, the FATHMM method has no 
obvious correlation with PhyloP, GERP++ RS, and LRT, while 
the correlation with MutationTaster is less significant. On the 
other hand, we observed significantly varied performances in 
detecting deleterious mutations by different methods (Table 4). 
Prediction errors by the in silicn methods highlight the limitations 
of these m(;thods and the need for cautious applications of the in 
silico prediction in data explanation. Among different methods, 
FATHMM has the lowest false negative rate (FNR = I %), but also 
the highest false positive rate (FPR = 95%) [Matthews correlation 
coefficient (MCC) = 0.127]. Considering the lack of correlation of 
FATHMM with the PhyloP, GERP++ RS, and LRT, caution 
should be taken when explaining the FATHMM results because of 
its high FPR and low MCC. The highest MCC scores were seen in 
the RadialSVM score (MCC = 0.474, FNR = 5%), PolyPhen-2 
HDFV (MMC = 0.447, FNR = 9%), PolyPhen-2 HVAR 
(MCC = 0.434, FNR =16%) and LR score (MCC = 0.393, 
FNR = 4%). 

Our investigation further disclosed significant differences of the 
quantitative performances of different methods, except SIFT, 
across different genes (Table 5). Varied performances across genes 
highlight another aspect of limitation of these in silico methods. 
The distribution of the prediction scores presented in Table 5 may 
be able to serve as a matrix to assist the assessment of functional 
effects of new mutations in these diabetes genes. 

The varied performances of these methods in different genes 
and the different scores of each method for different genes suggest 
that using gene-specific thresholds for deleterious mutations may 
improve the prediction performance of these in silico methods. We 
screened each gene and identified the gene-specific threshold with 
maximum MCC. Nonsynonymous single-nucleotide mutations in 
the diabetes genes from the NHLBI GO Exome Sequencing 
Project (ESP) [33,34], the ARIC samples [35] in the CHARGE 
Exome Sequencing Project [36], and the 1000 Genome Project 
[37], were used as controls without including mutations recorded 
in the HGMD database. Shown by our analysis (Table SI), we 
have been able to improve the prediction performance of each 
method in most cases, with the FATHMM method as an 
exception because of its nil/low FNRs in those diabetes genes. 
For example, the FNR of GERP++ RS for HNF4A mutations and 
the FNR of LRT for HNFIB mutations were decreased without 
any obvious change of their FPRs. On the other hand, redefined 
thresholds are able to decrease the FPRs of LRT for INS 
mutations, MutationTaster for ABCC8 mutations, LR score for 
INS mutations, LRT for ABCC8 mutations, MutationTaster for 
INS mutations, and MutationTaster for HNFIB mutations, 
without obviously increasing the FNRs. The general performances 
of different methods were summarized in Table 6. From low to 
high MCCs, the methods were sorted from left to right and from 
top to bottom. The average difference of MCCs and P value of 
each two methods was shown. 

The varied performance of different methods in difiFerent genes 
is related to specific molecular mechanisms of diabetes mutations. 
For the 41 INS mutations involved in this study, 34 mutations 
cause neonatal diabetes. These mutations exert diabetic effects by 
causing misfolding of the insulin protein, rather than inactivating 
the gene [49,50]. The dominandy inherited mode of the disease is 
from dominant negative mechanism, instead of haploinsufficiency. 
The misfolded insulin protein interferes cellular processes, leading 
to severe endoplasmic reticulum stress and potentially P cell death 



by apoptosis [50]. In contrast, a heterozygous individual with one 
copy of inactivating INS mutation may still have a sufficient 
response to metabolic regulation, thus without neonatal diabetes. 
For the prediction of neonatal diabetes mutations in the INS gene, 
a protein structure-based prediction method may thus have better 
performance than others. In this study, we see that PolyPhen-2 
with structure-based predictive features has better performance 
than the more sequence-based SIFT method (Table SI). Unlike 
other monogenic diabetes genes, the neonatal diabetes mutations 
in ABCC8 and KCNJll are gain-of-function mutations [51]. 
Sequence-based method like SIFT has also lower performance for 
these mutations than PolyPhen-2. 

We acknowledge the current publication bias of diabetes 
mutations (i.e. the bias towards identifying and reporting 
diabetes-causing mutations in the general human population). 
The diabetes mutations have been identified by studies involving 
much larger number of human individuals, while the genome 
sequencing projects involved limited number of human subjects. 
For a disease-causing mutation, no matter its low frequency, as 
long as the mutation is identified, it will be included. For example, 
in the case of GCK and HNFIA genes, the numbers of reported 
diabetes mutations are much larger than control mutations (479 
vs. 22, 324 vs. 78, respectively). We also want to emphasize the 
application of gene-specific mutations as functionally neutral 
controls. Our analysis showed that difiFerent methods using 
redefined thresholds by genome-wide control mutations, instead 
of gene-specific controls, tend to have poor performances (data 
available upon request). To acquire a satisfactory MCC tends to 
need a large number of both diabetes mutations and functional 
neutral mutations. The gene-specific prediction model proposed 
by our study will have further improved performance with the 
availability of sequencing data of a larger number of human 
individuals. 

In conclusion, the available in silico methods for the prediction 
of diabetes mutations have varied performances across difiFerent 
genes. In spite of the high statistical significance of correlations 
between difiFerent methods, difiFerent methods may not be able to 
replace each other. Because of variixi pc-rformances across genes, 
applying gene-specific thresholds when possible (i.e. for genes with 
a number of disease mutations identified and the ROC analysis 
feasible) may be able to increase the performance of prediction. 
For genes without sufficient numbers of mutations for the ROC 
analysis, a consensus threshold should be used [52]. Nevertheless, 
the limitations of the above methods warrant that new methods 
are being developed continuously. For example, Johansen et al. 
recently developed a sequence conservation-based artificial neural 
network predictor called NetDiseaseSNP [53]. Capriotti et al. 
developed a Meta-SNP algorithm for the detection of disease- 
associated nsSNVs, which integrates four difiFerent methods: 
PANTHER, PhD-SNP, SIFT and SNAP. They showed these 
methods are orthogonal with difiFerent biologically relevant 
relationships, and the integration of different methods achieved 
higher accuracy [54]. 

Supporting Information 

Table SI Method comparisons for gene-specific function 
prediction for non-synonymous mutations causing diabetes. 

(XLS) 
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